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Abstract 

Tomography is the three-dimensional reconstruction of an object from images taken at different 
angles. The term classical tomography is used, when the imaging beam travels in straight lines 
through the object. This assumption is valid for light with short wavelengths, for example in x- 
ray tomography. For classical tomography, a commonly used reconstruction method is the filtered 
back-projection algorithm which yields fast and stable object reconstructions. In the context of 
single-cell imaging, the back-projection algorithm has been used to investigate the cell structure or 
to quantify the refractive index distribution within single cells using light from the visible spectrum. 
Nevertheless, these approaches, commonly summarized as optical projection tomography, do not 
take into account diffraction. Diffraction tomography with the R.ytov approximation resolves this 
issue. The explicit incorporation of the wave nature of light results in an enhanced reconstruction 
of the object’s refractive index distribution. Here, we present a full literature review of diffraction 
tomography. We derive the theory starting from the wave equation and discuss its validity with the 
focus on applications for refractive index tomography. Furthermore, we derive the back-propagation 
algorithm, the diffraction-tomographic pendant to the back-projection algorithm, and describe its 
implementation in three dimensions. Finally, we showcase the application of the back-propagation 
algorithm to computer-generated scattering data. This review unifies the different notations in 
literature and gives a detailed description of the back-propagation algorithm, serving as a reliable 
basis for future work in the field of diffraction tomography. 
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Introduction 


Introduction 


Computerized tomography (CT) is a common tool to image three-dimensional (3D) objects like e.g. 
bone tissue in the human body. The 3D object is reconstructed from transmission images at different 
angles, i.e. from projections onto a two-dimensional (2D) detector plane. The assumption that the 
recorded data are in fact projections is only valid for non-scattering objects, i.e. when the imaging 
beam travels in straight lines through the sample. A common technique to reconstruct the 3D 
object from 2D projections is the filtered backprojection algorithm. The term classical tomography 
is commonly used to identify reconstruction techniques that are based on this straight-line approach, 
neglecting diffraction. 

When the wavelength of the imaging beam is large, i.e. comparable to the size of the imaged 
object, diffraction can not be neglected anymore and the assumptions of classical tomography are 
invalidated. The property that is responsible for diffraction at an object is its complex-valued 
refractive index distribution. The real part of the refractive index accounts for refraction of light 
and the imaginary part causes attenuation. Features in the refractive index distribution that are 
comparable in size to the wavelength of the light cause diffraction, i.e. the wavelike properties 
of light emerge. For example, organic tissue does not affect the directional propagation of x-ray 
radiation. Therefore, approximating the propagation of light in straight lines for CT is valid. How¬ 
ever, small refractive index changes within single cells lead to diffraction of visible light, requiring 
a more fundamental model of light propagation. Diffraction tomography takes wave propagation 
into account (hence the name backpropagation algorithm) and can thus be applied to tomographic 
data sets that were imaged with large wavelengths. Therefore, diffraction tomography is ideally 
suited to resolve the refractive index of sub-cellular structures in single cells using light from the 
visible spectrum. 

Because the refractive index distribution of single cells is mostly real-valued, they do not absorb 
significant amounts of light. Therefore, the main alteration of the wave front is due to refrac¬ 
tion which is revealed by a measurable phase change. Therefore, diffraction tomography requires 
imaging techniques that quantify these sample-induced phase changes, such as digital holographic 
microscopy (DHM). The general problem is depicted in figure 1. For different angular positions, 
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Figure 1, Tomographic data acquisition. An incident plane wave ito(r) is scattered by a transparent 
object with the refractive index distribution n{ r). A detector collects the scattered wave u(r). Multi-angular 
acquisition is facilitated by rotation of the sample. 
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1 Tomography without Diffraction 


images of single-cell sized objects are recorded at the detector. A plane wave uo(r) with a wave¬ 
length A in the visible regime propagates through a biological cell with a certain real refractive 
index distribution n(r). The recorded set of phase images, measured at different rotational posi¬ 
tions of the cell, is called a sinogram. Sinograms are the starting point for the refractive index 
reconstruction in three dimensions. 

The first section of the manuscript is a brief summary of non-diffraction tomographic methods. 
The two following sections introduce the wave equation and showcase the reconstruction from 
analytically computed sinograms. The subsequent chapters present the reconstruction algorithms 
[1-3] and their numerical implementation. The notation that we use here is similar to that used by 
the relevant literature (e.g. [4]). A list of symbols is given at the end of this manuscript. 

1. Tomography without Diffraction 

The first applications of computerized tomography (CT) used bone tissue as an inherent x-ray¬ 
absorbing marker. However, the marker can also be artificially introduced to the specimen. For 
example, in positron emission tomography (PET) radioactive tracers serve as markers for high 
metabolic activity. CT was also applied to biological specimens using wavelengths of the visible 
spectrum of light. The technique was termed optical projection tomography (OPT) [5]. In OPT, 
instead of measuring the absorption of the sample, the phase change introduced by the sample 
is measured. Thus, only the real part of the refractive index n(r) is reconstructed, whereas in 
classical x-ray tomography the imaginary (absorbing) part of the refractive index is measured. In 
OPT, fluorescent markers can be used to complement the measurement of refractive index, just 
as PET does for classical x-ray CT. The algorithms presented in this section do not take into 
account the wave nature of light. They are only valid in the limit of small wavelengths, i.e. x-ray 
radiation or for very small refractive index variations e n (r) = n(r) — n m of the sample n(r) from 
the surrounding medium n m . 

1.1. Radon Transform 

The Radon transform describes the forward tomographic process which is in general the projection 
of an n-dimensional function onto an (n — l)-dimensional plane. In the case of computerized 
tomography (CT), the forward process is the acquisition of two-dimensional (2D) projections from 
a three-dimensional (3D) volume 2 . The 3D Radon transform can be described as a series of 2D 
Radon transforms for adjacent slices of the 3D volume. For the sake of simplicity, we consider only 
the two-dimensional case in the following derivations. 

The value of one point in a projection is computed from the line integral through the detection 
volume [6]. The sample /(r) is rotated through cf >o along the y-axis. For each 2D slice of the 
sample /(r)| at y = y s , the one-dimensional projection p<p 0 (rn) = p^ 0 (xD,y s ) of this slice onto 


2 Keep in mind that in CT, the absorption of e.g. bone tissue is measured, whereas in optical tomography, the phase 
of the detected wave is measured. 
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1 Tomography without Diffraction 


a detector plane located at (xd, y s ) is described by the Radon transform operator R ^ 0 . 
P<fi 0 (xD,y s ) = ^0o{/(r)l y=ttl }(*D) 

= dvf(x(y), y v , z(v)) 


= J dv /(x D cos(f)o — v sin(/>o, y s: xp sin</>o + v cos 0o) (1.1) 

r xz = x +2 = x D + u ( 1 . 2 ) 

x'n = x cos 0o + z sin (po (1.3) 

v = — xsin(/»o + zcoscpo (1.4) 


Equation 1.2 defines a distance r xz from the rotational center of the 2D slice. Equations 1.3 and 1.4 
describe the dependency of the detector position at x = xd and the parameter v of the integral on 
the rotational position of the sample as defined by cj>o. The equations are illustrated in figure 1.1. 

Because the 3D Radon transform can be described by multiple 2D Radon transforms, the 3D 
reconstruction from a sinogram can also be described by multiple 2D reconstructions 3 . 



(a) 3D sketch 


(b) 2D slice integral pathway 


Figure 1.1, 3D Radon transform, a) Working principle of the three-dimensional (3D) Radon transform 
of a 3D object with the rotational axis y and the rotational angle </>q. For each slice of the object at y s (blue 
plane), a two-dimensional Radon transform is performed, b) The two-dimensional Radon transform at y s 
is computed by rotation of the object (white) through <p 0 (red coordinate system) and integration along v 
(green line) perpendicular to the detector line id- 


1.2. Fourier Slice Theorem 

The Fourier slice theorem is the central theorem in classical tomography. It connects the projec¬ 
tion data p^ 0 (x d) to the original image function f(x,z ) in Fourier space, which allows efficient 
reconstructions by means of the fast Fourier transform. 

li This is not valid for diffraction tomography as discussed in sections 4 and 5. 
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1 Tomography without Diffraction 


The projection of a two-dimensional image /(r) onto a detector line at an angle (fro can be written 
as the integral 


P<t> o(‘ t d) = Jdv f(x(v),z(v)). 


( 1 . 1 ) 


We define the unitary angular frequency Fourier transform of the one-dimensional data at the 
detector line with 


P(Uq (A'Dx ) 


1 


y/2n . 


'dx D Pcj >o (id) exp(— z^dx^d)- 


The actual Fourier transform F( k) of the two-dimensional image /(r) is given by 

1 


F(k x ,k z ) = 


2vr 


dxdz f(x, z) exp {—i{k x x + k z z)). 


(1.5) 


( 1 . 6 ) 


This formula can be rewritten as (subscript (fro denotes rotation; the Jacobian of a rotation is 1) 


F<f, 0 {k Dx ,k v ) = — Jjdxndvf < i )0 (xT>,v)exp(-i(kBxXn + Kv)) (1.7) 

Uo( x D,v) = f (xy) cos (fro - v sin (fro , xosin^o Tucos^o) ( 1 . 8 ) 

F^ (kr, x , K) = F(kr> x cos (fro - k v sin (fro, kn x sm(fro + k v cos(fr 0 ). (1.9) 

For the case k v = 0 (which implies slicing F( k) at the angle (fro), we arrive at [7-9] 

^. o (fcDx, 0 ) = —7==P,t, 0 {kn x )- ( 1 . 10 ) 

V Z7T 

This formula, known as the Fourier slice theorem, states that the Fourier transform of a projection 
at an angle (fro is distributed along a straight line at the same angle in the Fourier space of the 
image /(r). This theorem allows us to compute the inverse of the Radon transform operator by 
interpolating F(k) from P 0 O (/cdx) in Fourier space and subsequently performing an inverse Fourier 
transform. 

In order to compute the inverse Radon transform in Fourier space, we have to interpolate the 
data in Fourier space on a rectangular grid. However, this technique is afflicted with interpolation 
artifacts. Alternatively, one often uses the backprojection algorithm which is introduced in the 
next section. 


1.3. Backprojection Algorithm 

The backprojection algorithm connects the object function f(x,z) to the Fourier transform of the 
projection F 0 o (/cdx)- In order to derive it, we first express the object function f(x, z) as the inverse 
Fourier transform of F(k). 


f(x,z) 


— 11 dk x dk z F(k x , k z ) exp {i{k x x + k z z)) 


( 1 . 11 ) 
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1 Tomography without Diffraction 


We then perform a coordinate transform from (k x . kfi) to (A;d x > 4>o)- It can be easily shown that the 
Jacobian computes to 



k x — /tDx cos cj)Q ky sin cj )q 
k z = /cdx sin (j) 0 + ky cos o 
k v = 0 


( 1 . 12 ) 

(1.13) 

(1.14) 

(1.15) 


Therefore, combining equations 1.10 to 1.15, we get [7-11] 




1 

27T 



d<t >o | 7 'dx 


P<h (^Dx) 

V2tt 


exp[i/cDx(^ cos (j)Q + ^sin^>o)]- 


(1.16) 


Note that the integral over (f >o runs from 0 to tt. The integrals of k x , k z , and fcox are computed 
over the entire fc-space, i.e. over the interval (—oo, +oo). The term |A;dx| is a ramp filter in Fourier 
space which lead to the common term filtered backprojection algorithm in literature. However, we 
refer to it as the backprojection algorithm throughout this manuscript. 

Figure 1.2 depicts the sinogram acquisition and image reconstruction of a two-dimensional test 
target with the backprojection algorithm from 30 and 100 angular projections. Note that because 
of our chosen coordinate system, at the angle fio = 0, fco x coincides with the k x axis (xp <6 = ° x). 



(a) original image, (b) sinogram, 

500 x 500 pixels 500 projections 



(c) reconstruction (d) reconstruction 

from 30 projections from 100 projections 


Figure 1.2, Backprojection. a) The original two-dimensional image contains ellipses of different gray¬ 
scale levels, b) The sinogram shows 500 projection of image (a) from 0° to 180°. For the computation 
of the sinogram, only the circular region of the original image (red) was used, c) Reconstruction using 30 
equidistant projections, d) Reconstruction with 100 projections. The sinogram and reconstructions were 
created with the Python library radontea [12]. 


As discussed in the introduction, the backprojection algorithm assumes that light propagates 
along straight lines. It is thus not an optimal method for optical tomography which employs 
wavelengths that are in the visible spectrum of light. Thus, we need to include the wave nature of 
light in our reconstruction scheme. The next two chapters introduce the foundation of the Fourier 
diffraction theorem which yields better reconstruction results for objects that diffract light. 
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2. The Wave Equation 


2 The Wave Equation 


The propagation of light through objects in space r and time t follows the Maxwell equations. 
Light propagation through empty space can be described by the wave equation, which simplifies 
the vectorial description of electromagnetic fields to a scalar description of waves. When light 
propagates through inhomogeneous media, the wave equation does not anymore provide an exact 
description because the vectorial components of the electromagnetic field couple at e.g. refractive 
index boundaries. However, it is known that the wave equation is a good approximation for light 
propagation in inhomogeneous media [13] and we show that this scalar approximation of vectorial 
fields is valid for the discussed test targets. The wave equation reads 

■ va * M) - < 2 -« 

where T(r,t) is a time-dependent, scalar wave field in space, cq is the speed of light in vacuum, 
and n(r) is the spatial refractive index distribution. Note that is usually a constant coefficient 
which describes the speed of the propagating wave. Because we are not interested in the temporal 

information of T(r, f), we may simplify the wave equation by separation of variables. We then 

obtain the time-independent wave equation, known as the Helmholtz equation [14] 

(V 2 + k( r) 2 ) u(r) = 0, (2.2) 


where k{ r) is the wave number that depends on the local refractive index distribution n(r) and u( r) 
is the scattered field. To simplify the notation, we introduce the wave number inside the medium 
surrounding a sample k m and the local variation in refractive index inside a sample e n (r) as follows, 


n(r) 
k( r) 


27 rn m 

A 

n m + e„(r) 

nCr) 

A/ m 

Tl m 

k m (i + - 


n n 


(2.3) 

(2.4) 

(2.5) 


where A is the vacuum wavelength of the light. 


2.1. Homogeneous Helmholtz Equation 

For the homogeneous case, i.e. there is no sample (e^r) = 0), the wave equation becomes the 
Helmholtz equation for a homogeneous medium 

(V 2 + A; 2 ,) u Q ( r) = 0. (2.6) 

This second order ordinary differential equation has plane wave solutions of the form 

ii 0 (r) = a 0 exp(?X m s 0 • r), (2.7) 
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2 The Wave Equation 


where so is the normal unit vector and ao is the amplitude of the plane wave. Throughout this 
script, we use the convention that uo(r) in equation 2.7 defines a wave traveling from left to right, 
and thus uo(r) has a ‘+’-sign in the exponent. This convention is widely used in the literature 
dealing with diffraction tomography [1, 3, 4]. 


2.2. Inhomogeneous Helmholtz Equation 

Equation 2.2 can be rewritten as the inhomogeneous Helmholtz equation 

(V 2 + k 2 m ) u( r) = -/(r) u( r) 

'n(r) x 2 


with /(r) = k } 


n n 


- 1 


( 2 . 8 ) 

(2.9) 


In order to deal with the inhomogeneity /(r), also called scattering potential, we can make use 
of the Green’s function G(r — r') that is defined by equation 2.10 4 . The Green’s function for the 
homogeneous Helmholtz equation is shown in equation 2.11°. 

(V 2 + k^) G( r - r') = —<5(r - r') 

exp(ifc m |r — r'|) 


G( r - r) = , , 

47t | r — r' | 

Here, 5(r — r') is the Dirac delta function with the translational property 


( 2 . 10 ) 

( 2 . 11 ) 


Jd 3 r’ 8{r - r') g( r') = g( r). 


( 2 . 12 ) 


By knowledge of the Green’s function we can derive an integral representation for the scattered 
held u(r). Using equation 2.12, the right side of equation 2.8 can be integrated to include the 
Green’s function (eq. 2.10) 6 


/(r) it(r) = j d 3 r 5(r — r) f(r) u(r) 

= ~ jd 3 r' (V 2 + km) G(r - r') f(r') u(r') 
= — (V 2 + k^) jd 3 r G{r — r') /(r') u(r'). 


(2.13) 

(2.14) 

(2.15) 


The comparison of equations 2.8 and 2.15 suggests that the scattered held u( r) can be written as 

u(r) = /dVG(r-rWMr') (2-16) 

which is the integral equation to the inhomogeneous wave equation (eq. 2.8). 

4 The Green’s function is defined as the solution to the inhomogeneous problem: When the operator (V 2 + k is 
applied to G(r — r'), one obtains the Dirac delta distribution <5(r — r'). 

5 Note that our notation implies the (+)-sign and a normalization factor of 47r for the Green’s function. For a 
derivation see for example [15], section 2.4 or [16], section 7.2. 

6 Note that the Green’s function depends on the absolute value of |r — r'| and thus 
(V' 2 + kl) G(r - r') = (V 2 + k 2 m ) G (r - r') 


The Theory of Diffraction Tomography 


9 



3 Approximations to the Scattered Wave 

3. Approximations to the Scattered Wave 

Equation 2.16 has no analytical solution. However, under certain conditions approximations can 
be applied to find a solution. In this section, we derive the Born and Rytov approximations. Both 
approximations define the scattered field u(r) as a superposition of the incident plane wave «o( r ) 
and a scattered component u s ( r). 

u(r) = u 0 (r) + u s (r) (3.1) 

3.1. Born Approximation 

The Born approximation uses the property of the homogeneous component 

(V 2 + /4) uo(r) = 0 (2.6) 

to rewrite the inhomogeneous wave equation for the scattered component u s (r) 

(V 2 + A&) u(r) = (V 2 + O 'Ug(r) = -/(r)u(r). (3.2) 

By once again using the translational property of the delta function (eq. 2.12) we obtain an iterative 
equation for u s ( r) 

u s ( r) = Jd 3 r' G(r — r') f(r r ) u(r') (3.3) 

and we can write the Lippmann-Schwinger equation for the field u{ r) with the perturbation /(r) 
[14] 

u(r) = tto(r) + Jd 3 r' G(r — r') f(r') u(r'). (3.4) 

By iteratively replacing u( r) in the above integral, one obtains the Born series. The Born approx¬ 
imation is actually the first iteration of the Born series [14]. 

. x Born / \ / \ / \ 

u[y) « ^ 0 (r)+^B(r) (3.5) 

ub( r) = Jd 3 r' G( r - r) /(r') u 0 (r'). (3.6) 


Validity of the Born Approximation 

In the derivation above, approximating the scattered component u s ( r) as the first term of the Born 
series «B(r) implies that 

Jd 3 r' G( r - r) /( r) [rto(r') + u s ( r')] « jd 3 r' G( r - r) /( r') u 0 (r') (3.7) 

or u 0 (r) +u s (r) « tt 0 (r). (3.8) 

Thus, the Born approximation holds for the case U u s (r) <C ito(r)”, i.e. the contributions of ampli¬ 
tude and phase of the scattering component u s ( r) are small when compared to the incident plane 
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3 Approximations to the Scattered Wave 


wave «o( r )- Since both it s (r) and uo(r) are complex-valued functions, the above relation implies 
that (i) there is only little absorption by the specimen and that (ii) the overall phase change AT 
must be small. Because we are interested in biological cells whose imaginary part of the refractive 
index is insignificant, (i) absorption is negligible and we thus focus on (ii) the phase change AT in¬ 
troduced by the cell. The phase change AT that a scattered wave u s (r) experiences when traveling 
through a sample along a certain path parametrized by s pa th can be approximated with ray optics 


A _ 2tt 
AT « — 
A 


J rfSpath n( r) - StotUm , 

\^path 


(3.9) 


where n(r) is the refractive index distribution inside the sample, s is the approximate thickness 
of the sample, and n m is the refractive index of the medium surrounding the sample. For a 
homogeneous sample with the refractive index n s , the absolute phase change computes to 


A , 27T 27T 

AT = — s(n s - n m ) = —se n 


(3.10) 


This equation can be interpreted as a comparison of the phase change AT over a period of 27 t with 
the change of the optical path length s(n s — n m ) over one wavelength A. 

We now want to write equation 3.10 in its differential form. We express the total differential of 
AT with respect to the spatial distance s and the refractive index variation e n as 


d( AT) e n s 

-*r= a* + a*' 


(3.11) 


The phase change AT can have two contributions, namely the thickness of the sample 


d( AT a ) _ 6n 
2 tt “A 

and the refractive index variation inside the sample 

d(AT B ) 

2t r “ A 


(3.12) 


(3.13) 


Note that in general e n is dependent on r (3D) and that AT is measured at the detector plane 
td (2D). If the refractive index of the sample is fixed, then the total phase change solely depends 
on the thickness of the sample. If the thickness of the sample is fixed, the local variations inside 
the sample determine the absolute phase change. If we wanted to consider a phase change that is 
introduced by a refractive index variation over the distance of one wavelength A, then we would 
have to set s = A. 

With these considerations, we can answer the question of the validity of the Born approximation. 
We interpret the inequality “tt s (r) <C ito(r)” as a general restriction: The overall phase change 
must be much smaller than 27T. According to equation 3.10, this translates to a restriction for the 
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3 Approximations to the Scattered Wave 

optical path difference s(n s — n m ) which must be much smaller than the wavelength A of the used 
light 


A<f* < 2vr (3.14) 

s(n s — n m ) <C A. (3.15) 

The fact that the Born approximation only allows to observe optically thin samples is a serious 
drawback. A second approach to the problem is the Rytov approximation. 


3.2. Rytov Approximation 

In order to derive the Rytov approximation we assume that the scattered wave u( r) and the incident 
wave rto(r) have the form 


u{ r) = exp(</?(r)) = u 0 ( r) + u s { r) (3.16) 

u 0 (r) = exp(y? 0 (r)) (3.17) 

T’(r) = <£ 0 (r) + ^s(r) (3.18) 

where we use the complex phases 

<^(r) = i<3?(r) + ln(a(r)) (3.19) 

T’o(r) = *^o(r) + ln(a 0 (r)) (3.20) 


to denote the complex exponent containing phase and amplitude of the wave function. Note that 
u s ( r) is now computed from y? s (r) with 


u s { r) = u(r) - u 0 (r) 

= exp(^ 0 (r)) [exp(y? s (r)) - 1] . 

Using the equations above, the inhomogeneous wave equation becomes 

(V 2 + /4M r) = —/(r)u(r) 

(V 2 + fc 2 J exp(<^(r)) = -/(r) exp(y>(r)). 

We can then compute the term V 2 exp(<£>(r)) 

V 2 exp(y?(r)) = V [exp(</?(r)) • Vy?(r)] 

V 2 exp(<£>(r)) = exp(</?(r)) vV(r) + (V</?(r)) 2 

to obtain a differential equation for g?(r). 

exp(v?(r)) V 2 </?(r) + (V^r)) 2 + A& = -/(r) exp(<^(r)) 

VV(r) + (V#)) 2 + e = -/(r) 


(3.21) 

(3.22) 


(3.23) 

(3.24) 


(3.25) 

(3.26) 


(3.27) 

(3.28) 
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3 Approximations to the Scattered Wave 

Equation 3.28 is a non-linear differential equation for the complex phase <£>(r). In the same manner, 
a differential equation for <^o( r ) can be derived 

VVo(r) + (VVo(r)) 2 + *4 = 0. (3.29) 

The next step is to insert equation 3.18 into equation 3.28 to find a differential equation for <y? s (r). 

V 2 [<A)(r) + ¥>s(r)] + (V[<^ 0 (r) + ^ s (r)]) 2 +*4 = -/(r) (3.30) 

- '-v-' — 

(V</>o(r)) 2 +2Vy?o(r)-V¥Js(r)+(Vv3s(r)) 2 


The terms marked with an underline compute to zero (eq. 3.29) and the equation above becomes 

VV s (r) + 2Vy? s (r) • Vv? 0 (r) + (V(^ s (r)) 2 = -/(r). (3.31) 

To simplify this expression we need to consider 

V 2 u 0 (r)<^ s (r) = V 2 it 0 (r) -^ s (r) + 2 Vit 0 (r) -V^ s (r) + ii 0 (r)VV s (r). 


- fc m«o(r) 


u 0 (r)Vv5o(r) 


I 


(3.32) 


(3.33) 


(V 2 + /4)u 0 (rV s (r) = 2u 0 (r)V</? 0 (r) • Vy? s (r) + u 0 (r)V 2 (/? s (r) 

If we multiply equation 3.31 by uo(r) then we can substitute with equation 3.33 to obtain 

(V 2 + *4)«o( r ) <p B (r) = -u 0 (r) [(V<^> s (r)) 2 + /(r)]. (3.34) 


Rytov 




Rytov 


/( r) 


The Rytov approximation assumes that the phase gradient V<^ s (r) is small compared to the per¬ 
turbation /(r). We can now make use of the Green’s function again (eqns. 2.10 to 2.16) and arrive 
at the formula for the Rytov Phase Vhi( r ) [4]. 


u 0 (r)^ R (r) = Jd 3 r' G( r - r) /(r) u 0 ( r) 
Jd 3 r' G (r — r') f(r') uq( r') 


<pr(*) = 


uo(r) 


(3.35) 

(3.36) 
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By comparing this expression to the first Born approximation (eq. 3.6) we find that we can compute 
the Rytov approximation rt R (r) from the Born approximation «B( r ) and vice versa. 


7ht(r) 

_ UB{r) 

«o(r) 


(3.37) 

HR(r) 

= «o(r) 

[expf nB(r)N ) ll 

r p v«(r)i 

(3.38) 

«B(r) 

= “° (r)ln (-o(r) +1 ) 

(3.39) 


= «o(r)^R(r) 

(3.40) 

u( r) 

Born 

') + u B {r ) 


~ ^o(r 


u(r) 

Rytov 

~ n 0 ( 

r) + « R (r) 



Note that the complex Rytov phase </? R (r) also contains the amplitude information of the scattered 
wave. In section 5.5 we introduce a backpropagation algorithm for the Born approximation. This 
algorithm may also be used in combination with the Rytov approximation by using equation 3.39. 
In practice, calculating the logarithm of exp( 9 ?(r) — <^o( r )) involves calculating the logarithm of the 
amplitudes and the phases. 

ln(exp(<^(r) - <^ 0 (r))) = In + i ($( r ) “ $o(r)) (3-41) 

Because the phase <h(r) — <f>o(r) is modulo 27r, a computational implementation of the Rytov 
approximation must contain a phase-unwrapping algorithm [17]. 

Validity of the Rytov Approximation 

From our approximation in equation 3.34, we can tell that the Rytov approximation is valid for 
the case 


(V^ s (r )) 2 < /(r) 

eq. 2.9 

< k 


2 

m 



n(r ) 2 > n 2 m 


~ (VPs(r )) 2 
1.2 


(3.42) 

(3.43) 


We are interested in a validity condition that connects the refractive index with its gradient. We 
insert the definitions of the wave vector k m and the refractive index distribution n(r) (eq. 2.3, 2.4) 
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to retrieve a condition for the variation in refractive index e n (r). 


2 ^ 2 /V 9 ? s (r)A 


n(r) J > n 2 


2im n 


+ n n 


n(r) 2 - r? m > 
e n (r) 2 +2n m e n (r) > 


V^s(r)A 

2vr 

V^s(r)A 

2vr 


(3.44) 

(3.45) 

(3.46) 


Because the local variation e n (r) is small, we may neglect' e n (r) 2 . The resulting constraint for the 
phase gradient is [4] 

IV</? s (r)| ^ yj2n m |e n (r)| 

2vr ^ A 
|dv? s (r)| ^ yj2n m |e n (r)| • |dr 
2vr ^ A 


(3.47) 

(3.48) 


For any position r inside a sample, equation 3.48 reads: 

The sample induces a phase change over a period of 2ir radians. This number must be smaller 
than the variation in refractive index e n (r) along the corresponding optical path scaled by the used 
wavelength A. 

Note that the total phase gradient is computed from 


|V</?(r)| = |\Vo(r) + V</? s (r)| (3.49) 

|V</?(r)| = |ik m + VV s (r)|. (3.50) 


Thus, compared to the Born approximation, where the overall phase change must be smaller than 
27 t, the Rytov approximation is also valid for optically thicker samples. 

The Rytov approximation is accurate for small wavelengths and breaks down for large variations 
in the refractive index [18]. The following calculations attempt to derive a statement of validity for 
the Rytov approximation that only depends on the refractive index, which is difficult considering 
its non-linear, but exponential description of wave propagation. When we insert equation 3.13 into 
equation 3.48 (A4>B(r) = ^ s ( r ))i we obtain the restriction for the Rytov approximation. 


de n (r) y / 2re m |e n (r)| 

| dr | s 


(3.51) 


Here, we assumed that we can replace a change in the local complex phase dtp(r) with a variation in 
the local refractive index de n ( r). This estimate is valid when light propagates approximately along 
straight lines, as described by equation 3.13. Furthermore, this estimate does not cover scattering 


7. 


This can be shown by solving the quadratic equation 3.46 for e n (r) and Taylor-expanding for small V<p s (r) to the 
second order. 
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4 Two-Dimensional Diffraction Tomography 


at small objects and thus, s > A. The quantity s can be interpreted as a characteristic length scale 
within the imaged object below which light propagation can be approximated along straight lines 
(eq. 3.13). 


|Vn(r)| <C 


y / 2 n m |ra(r) - 


s > A 


(3.52) 


The validity of the Rytov approximation is not dependent on the absolute phase change introduced 
by the sample, but on the gradient of the refractive index within the sample. This makes the Rytov 
approximation applicable to biological cells. The application of the Rytov approximation to the 
three-dimensional backpropagation algorithm is described in secion 5.3. 


4. Two-Dimensional Diffraction Tomography 

The name Fourier diffraction theorem was introduced by Slaney and Kak [4] in the 1980’s. In the 
limit of small wavelengths, the Fourier diffraction theorem in the Rytov approximation converges 
to the Fourier slice theorem (section 5.3, [2]). In this section, we derive the two-dimensional 
backpropagation algorithm, as described by Kak and Slaney [4], which was implemented in the 
C programming language in 1988 s . At the end of the section, we test the backpropagation algorithm 
with artificial data that were computed using Mie theory. Mie theory provides solutions to the 
Maxwell equations in the form of infinite series for scattering objects consisting of simple geometric 
shapes. It is thus ideally suited to test the approximations made in sections 2 and 3. The 3D 
theory can be derived analogous to the 2D theory. It is described and discussed in greater detail 


in section 5. Throughout this section we use two-dimensional vectors, i.e. r = (x, z ). 

4.1. Fourier Diffraction Theorem 

We start from the inhomogeneous wave equation as discussed in section 2.2. 

(V 2 + A&) u( r) = -f(r)u(r) (4.1) 

The Green’s function in the two-dimensional case is the zero-order Hankel function of the first kind. 

(V 2 + kl) G(r - r') = —<5(r - r') (4.2) 

G(r - r') = l -H^\k m |r — r'|) (4.3) 

H^\k m |r - r'|) = ^ Jdk x - exp{i [k x (x - x') + 7 (z - z')] } (4.4) 

7 = vT 2 i - kl (4.5) 


We observe restrictions for the Cartesian coordinates for the vector so describing the incoming 
plane wave and the vector s describing an arbitrarily scattered wave. The following substitutions 


implemented by Malcolm Slaney, available at http://slaney.org; see also [19] 
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are defined by our notation: 


/c x — k m p , 

7 = k m M 

(4.6) 

M = y/l - p 2 , 

M 0 = \Jl ~ Po 

(4.7) 

s = (p, M ), 

So = (po, Mo) 

(4.8) 


The parameters p and po describe the fc x -component of the vectors s and Sq. They must fulfill the 
relation 0 > p,po > 1. In order to have the data consistent with section 1 (Fourier slice theorem), 
so points into the ^-direction when <f >o = 0. This means that so = (— sin</>o, cos^o)- We rewrite 
the Green’s function accordingly. 

G( r - r) = ^- Jdp exp {ik m [p(x - x) + M(z - z')\ } (4.9) 

In our notation, the incoming plane wave uo(r) with amplitude ao, propagation direction So, and 
wavenumber k m can be written as 


wo(r) = a 0 exp(?X’ m s 0 r). 

The first Born approximation in two dimensions then reads (see section 3.1). 

u B (r) = [ fd 2 r'G(r - r')f(r')u 0 (r') 


We define the two-dimensional Fourier transform F(k ) of the scattering potential /(r) 

F( k) = JJd 2 r /(r) exp(—ikr) 

/(r) = JJd 2 kF(k)exp{ikr) 

and we keep in mind the identity of the Dirac delta distribution for any given p and a 

S(p — a) = — [dx exp (i(p — a)x). 

2vr J 

We now insert equation 4.9 into equation 4.11. 

“ B(r) = i Ih' J dp S e ^ ik ^ x - x,)+ - Z,) F 

x/(r / )a 0 exp(ifc m (p 0 ^ / + M 0 z ')) 

The integral over r' can be replaced with the Fourier transform of 

F(k m (s - s 0 )) = If d 2 r'f(r) exp(-?X m (s - s 0 )r) 


(4.10) 

(4.11) 

(4.12) 

(4.13) 

(4.14) 

(4.15) 

(4.16) 
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4 Two-Dimensional Diffraction Tomography 


The equation for the Born approximation now reads 


/ \ m /* 2tt ~~. j , ,. . i 

«B(r) = — dp —F(k m (s — s 0 )j exp(?A: m sr). 


(4.17) 


We measure the field at the detector line r —> ro = (*dPd) and at the angle (fro- Here, Id is the 
distance of the detector plane from the rotational center of the sample. 


«B,</.o( r D) = — dp —F(k m (s - s 0 )) exp(ik m srr,) 


(4.18) 


The subscript (fro denotes the angular position of the sample and the direction of the incoming 
plane wave with respect to the sample. In order to arrive at the Fourier diffraction theorem in two 
dimensions, we perform a one-dimensional Fourier transform of its (ro) along xd- 


tWo^’Dx) = Jdx D Jdp^F(k m(s So)) exp(?X’ m (px D T T7/p>))x 

x exp(-ifc Dx x D ) 

We now identify the delta distribution 


5(k m p - fc Dx ) = Jdx D exp (i(k m p - fe Dx )a: d) 


which in turn we can use to solve the integral over dp. 

Note that 5(k m p - fc Dx ) = ^ S(p - k Dx /k m ). 

UB,<f> 0 {kD x ) = J d P - s 0 )) exp(ik m MlD)6(k m p - A: Dx 

Solving the integral implies replacing all occurrences of p with . 


^B,0o(^T> x ) = 


\/2irk n 


-F(k m (s - s 0 )) exp ik m \ 1 - 




Finally, we arrive at the 2D Fourier diffraction theorem 


^B,^o(^Dx) = 


iao n r 1 


fc m V 2 M 


-F(k m (s - s 0 )) exp(ik m Ml B ) 


(4.19) 


(4.20) 


(4.21) 


(4.22) 


(4.23) 


Note that we have used the substitution M = \ 1 — (to simplify the expression. Solving for 


the Fourier transformed object F yields 


F(k m (s - s 0 )) = - \ -M£/b^ 0 (£;d x ) exp 

V vr ao 


(4.24) 
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The restriction in equation 4.7, rewritten in equation 4.26, forces the one-dimensional Fourier 
transform of the scattered wave £/b ,</> 0 (&Dx) to be placed on circular arcs in Fourier space: 


&: m s = (fcox cos (f >o — k m M sin (f> o, /cd x sin (j)Q + k m M cos (j>o) 
km M = \Jk £ - kl~ 


(4.25) 

(4.26) 


The argument k m ( s — so) shifts the circular arcs in Fourier space such that ?7b,</> 0 (0) is centered at 

f(o,o). 


Comparison to the Fourier Slice Theorem 

Equation 4.24 describes the Fourier diffraction theorem in 2D. When compared to the Fourier 
slice theorem from equation 1.10 in section 1.2, a few differences become apparent. We can write 
equation 4.24 in the same manner as equation 1.10, with the subscript <j )o denoting the rotation of 
the object /(r). The main differences are a complex factor and a different distribution of the data 
in Fourier space, as illustrated in table 4.1. 
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Fj, 0 ( 5 ) — -A 

• o(^Dx) 


Fourier Slice Theorem 

(equation 1.10) 

Fourier Diffraction Theorem 

(equation 4.24) 

Sinogram 

P (&Dx) 

Fourier transform 
of projections (fc Dx ) 

Fourier transform of complex 
scattered waves Cb,^ 0 (&Dx) 

Factor 

A 

A = 1 

A =- m exp(— ik m Ml^>) 

do 

Coordinates 
(fc x , kz) 
sliced at </> o 

k x = kn x 

k z — k t — 0 
(straight line) 

k x = kn x 

k z = \J k^ - k^ x - k m 

(semicircular arc) 


Fourier 
space F(k) 
coverage 
(180°) 




Table 4.1, Classical versus diffraction tomography. The equation above the table combines the Fourier 
slice theorem and the Fourier diffraction theorem, pointing out their common structure. However, there are 
two major differences between the theorems. (1) Diffraction tomography data are multiplied by a complex 
factor A/l. (2) For diffraction tomography, the Fourier transform of the obtained data is distributed along 
circular arcs (not straight lines) in Fourier space. The plots show a visualization in Fourier space for data 
that are sampled at a frequency of 1 /k la for angles between 0° and 180°. 


The Theory of Diffraction Tomography 


20 


















4 Two-Dimensional Diffraction Tomography 

4.2. Backpropagation Algorithm 

The term backpropagation comes from an interpretation of the mathematical formalism which is 
similar to the backprojection formula derived in section 1.3. The backpropagation algorithm [4] is 
a solution to the inverse scattering problem. The reconstructed scattering potential /(r) is directly 
computed from the input data Ub (k). For this, we need to perform a coordinate transform from 
(k x ,k z ) to ( k]Oxi )• We start by computing the inverse two-dimensional Fourier transform of 
equation 4.24. 


^ / 2 ik ^ 

F(k m (s - s 0 )) = - \ - -MUb,^ 0 (A:dx) exp(-f£: m Af7 D ) (4.24) 

V 7T ao 

/(r) = ~~ J jdk x dk z MF B ,^ 0 (fex) exp(-ik m Ml B ) x 

x exp(i/c m (s — so)r) (4.27) 

(k x , k z ) =k m {s - s 0 ) (4.28) 


As described in table 4.1, the input data are distributed along circular arcs in Fourier space. The 
orientation of these arcs is defined by the acquisition angle 0 o and is described by D ^ 0 . 


n (COS 00 

^ sin 0 o 

k = Dfo k' 


— sin 0 o 
cos 00 


(4.29) 

(4.30) 


Here, k denotes the non-rotated Fourier space, whereas k' denotes the positions of acquisition at a 
certain angle 0o- At 0o = 0, k and k' coincide. We have defined the angle 0o such that, k x = fco x 

and therefore k' z = \Jk^ — k^ x — k m . The coordinate transform from (k x ,k z ) to (&Dx,0o) is now 
fully described. 


k z 


fcox COS 00 



fc Dx sin 00 + 



sin 0 o 

(4.31) 

COS 00 

(4.32) 


In order to replace the integrals over k x and k z with the integrals over k^x and 0o, we compute the 
Jacobian matrix J and its determinant. 


dk x dk z 
dkn x d<j) o 


(cos 00 + , 2 P _ X 3 sin 00 

V k Dx 

^sin 0 o - cos 0 o 


-fc Dx sin 0 o 




cos 00 


k~Dx COS 00 




sin 0 o 


(4.33) 

(4.34) 
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The determinant of the Jacobian J computes to 


det(J) = fc Dx 



^m^Dx 

\J k m ~ ^Dx 


^’m^Dx 



(4.35) 

(4.36) 


We insert the coordinate transform into equation 4.27 and obtain the backpropagation formula. 


/ ( r ) 




MUB,<t> 0 (kr> x ) exp(-ik m Ml D ) x 


x exp(ifc m (s — so)r) 


(4.37) 


Note that the integration over 0o goes from 0 to 27 t. This is necessary because the Fourier space 
needs to be covered homogeneously by the projection data. A coverage of only 180° leads to an 
incomplete coverage in Fourier space as depicted in table 4.1. This aspect is important to keep in 
mind for later experimental realization. The necessary double-coverage in Fourier space leads to 
the additional factor 1 / 2 . Furthermore, we express (s — so) in terms of a lateral (tj_) and an axial 
(so) component (eq. 4.31 and 4.32). 

k m (s - s 0 ) = k Dx t ± + k m (M - 1) s 0 (4.38) 

s 0 = ( po , Mo) = (— sin0o, cos0 o ) (4.39) 

t± = (—Mq, po) = (cos0o, sin0 o ) (4.40) 


! 

By assuming that ( k m M ) 2 = k 2 n — k > 0, we can rewrite the backpropagation formula as 


/( r) = 


ik m 

ao(2n) 3 / 2 



r2ir 


#o I^Dxl t/B, 0 o(^Dx)exp(-iA: m M/ D )x 

x exp[?'(fc Dx tj. + k m (M - 1) s 0 )r] . 


(4.41) 


Figure 4.1 illustrates the acquisition and reconstruction process for a non-centered cylinder. The 
data were computed according to Mie theory. Note that the reconstruction with the Born approx¬ 
imation is not resembling the original structure of the cylinder and that the reconstruction quality 
with the Rytov approximation is dependent on the total number of projections that were acquired 
for the sinogram. Also note that the main contribution to the reconstructed refractive index dis¬ 
tribution comes from the measured phase. The measured amplitude only slightly influences the 
reconstruction. Details of the simulation and reconstruction process can be found in figure 4.2. 
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(d) Born (250 projections) (e) Rytov (50 projections) (f) Rytov (250 projections) 

Figure 4.1, 2D Backpropagation. Image (a) shows the refractive index distribution of a dielectric 
cylinder. White indicates the refractive index of the medium and black the higher refractive index of the 
cylinder. The forward-scattering process was computed according to Mie theory. Because the cylinder is not 
centered, the characteristic shape of the sinogram is visible in amplitude (b) and phase (c). The lower images 
show the reconstruction of the refractive index with (d) the Born approximation using 250 projections, (e) the 
Rytov approximation using 50 projections, and (f) the Rytov approximation with a total of 250 projections. 
See figure 4.2 for details. 
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cross-section [A] 

Figure 4.2, 2D Backpropagation cross-sections. The refractive index of the medium is n m = 1.333 
and the local variation inside the cylinder is e„(r) = n(r) — n m = 0.006. The radius of the cylinder is 30A 
(vacuum wavelength A). The scattered wave is computed according to Mie theory at an optical distance of 
Id = 60A from the center of the cylinder and sampled at A/2 over 250 pixels. The Mie theory computations 
are based on [20]. The refractive index map is reconstructed on a grid of 250 x 250 pixels from 250 projections 
over 27 t. The reconstruction was performed with the Python library ODTbrain [21]. 
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Comparison to Backprojection 

The main differences between the backprojection and backpropagation algorithms (dependencies 
on Id and sq) are summarized in table 4.2. 



f{x ’ V) =(2*)^I dkD 'L 

271 ^ 

#o exp (iB) |/c Dx | P<t> o(fcDx) 


Backprojection 

(equation 1.16) 

Backpropagation 

(equation 4.41) 

Sinogram 

Pft o (^Dx) 

Fourier transform of projections 

P (^Dx) 

Fourier transform 
of complex scattered waves 

^B,</>o(^Dx) 

Factor 

A 

A - I 

A 2 

(double coverage) 

i/c m 

ao 

Exponent 

B 

B = fc Dx (tj_r) 
t± = (cos 0 o , sin 0 o) 

B = —k m MlD + &Dxt-i_r + 

+ k m (M - l)s 0 r 
t_L = (cos 0o, sin 0 O ) 
s 0 = (-sin 00 , cos 0 o) 


Table 4.2, Backprojection versus backpropagation. The equation above the table illustrates the 
similar structure of the backprojection and the backpropagation formula. Nevertheless, the backpropagation 
formula contains an exponent B , resulting from the first Born approximation, that depends on Id and So 
which increases the complexity (See also table 4.1). Note that the backprojection formula has a factor of 
!/2 due to the double coverage over (j) 0 . The necessity of the 27r-coverage (360°) for the backpropagation 
algorithm is illustrated by the visualization of the Fourier slice and diffraction theorems in the figures of 
table 4.1. 


5. Three-Dimensional Diffraction Tomography 

Two-dimensional (2D) diffraction tomography is valid in three dimensions (3D) only for infinitely 
elongated objects (i.e. cylinders). For objects that exhibit inhomogeneities along the third dimen¬ 
sion, a 3D theory needs to be applied. The 2D and the 3D theory are very similar. There are 
only a few dimension-specific differences like the position of the Fourier transformed fields along 
spherical surfaces (3D) instead of circular arcs (2D). This section uses the previous sections as a 
basis to introduce and discuss the 3D version of diffraction tomography. The notation stays the 
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same, except that a third axis y is introduced in every vector. As in section 4, we conclude by 
testing the 3D backpropagation algorithm with data computed according to Mie theory. 

5.1. Fourier Diffraction Theorem 

There are at least two ways to derive the 3D Fourier diffraction theorem. The first uses a double 
integral representation of the Green’s function. The second makes use of the convolution theorem 
that connects the convolution of two functions with their product in the Fourier domain. We will 
thoroughly investigate the first and only briefly discuss the second method. 

The Green’s function (equation 2.11) of the Helmholtz equation can be rewritten using the double 
integral representation as shown by Banos et. al 9 ([22], section 2.11). 

G(r - r') = ^ JJ dpdq exp {ik m [p(x - x') + q(y - y') + M(z - z ')] } (5.1) 

M = y/l -p 2 - q 2 (5.2) 

s =(p,q,M) (5.3) 

Here, we define £: m s as the wave vector of a plane wave with the wave number k m . Note that we 
have introduced the coordinate q in addition to p from the 2D case. In order to keep M real, we 
now have the restriction p 2 + q 2 < 1. If this restriction was violated, we would allow evanescent 
waves which would complicate the inversion process [4], An incoming plane wave has the normal 
vector sq that is defined by 


it 0 (r) = a 0 exp(ifc m s Q r) (5.4) 

with s 0 = (p 0 ,q 0 ,M 0 ). (5.5) 

In section 3.1 we showed that the Born approximation of u s (r ) reads 

u B (r) = JJJd 3 ' 1 '' G ( r ~ r ')f( r ') u o(r')- (5.6) 

We define the unitary Fourier transform 10 F’(k) of /(r) in 3D as 

^( k ) = ^ Jfjd 3 r /(r) exp(-i kr) (5.7) 

= ( 2 k ) 3 / 2 J[J d3k 6Xp ^ kr )‘ ^ 5 ' 8 ^ 


Furthermore, we define the position of the detector such that its normal vector np is parallel to the 
incident plane wave vector k m so (np = so). This is a valid assumption, because we only rotate the 
cell and thus the spatial arrangement of incoming plane wave uo(r) and detector do not change. 

u Banos’ notation does not include the prefactor l/(47r) in the Green’s function (see eq 2.11). Therefore, the prefactor 
in equation 5.1 is l/(87r 2 ) and not l/(27r). 

1(, Please note that other authors, e.g. Devaney [2] may use the non-unitary Fourier transform. The prefactors 
(l/(2-7r) 3//2 ) may differ from our notation. 
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In the derivations that follow, we make use of the following definition of the Dirac delta function. 


5(x — a) = — [dp exp (ip(x — a)) 
2vr J 


(5.9) 


We insert the definitions of the Green’s function and the incident plane wave into the formula for the Born approximation 
M r ) = Z g^ 2 m JJJ d 3 r' j J dpdq exp {-ik m [(p - p 0 )x' + (q - q Q )y' + (M - M 0 )z'] } ■ exp (ik m (px + qy + Mz) ). 

sr 

(5.10) 


to 

-4 


We may interpret the integral over r' as the Fourier transform of /(r) that is shifted by —k m sq in Fourier space 


M r ) = l ~~^T ( 27r ) 3/2 JJ dpdq -jjF( k m (p - p 0 ), k m (q - go), k m (M - M 0 ) ) • exp (ik m (px + qy + Mz)). (5.11) 

fcm(s-So) 


Any vector ro = (xb,Vt>, zt>) in the detector plane is defined by np(r — ro) = 0, where np is the normal vector on 
the detector plane. We previously placed our detector according to no = So, which implies a rotation of the sample 
perpendicular to the propagation direction so of the incident plane wave ito(r). Thus, z d = Id is a constant and describes 
the distance of the detector from the center of the rotation axis. The detected field at the detector plane is thus 


M r ) Id = 


ITT CIq kiYi 

( 27 r ) 3 / 2 


dpdq —F(k m (s - s 0 )) • exp (ik m (px-D + qy d + Ml D )). 


(5.12) 


Now we compute the two-dimensional Fourier transform of the recorded data UB, 0 o ( r D) in the detector plane (xd, 2 /d)- 


^B,0 o ( k D) = 


ITT Uq k^i 
(2tt) 5 /2 


dx-odyn / /dpdq —F(k m (s - s 0 )) • exp (ik m (pxv + qy D + Ml D )) • exp(-i(fc Dx x D + &DyM 


(5.13) 


Note that since k^ x + k^ y + k^ z = kf n . i.e. there is no inelastic scattering, ( k D) must be on a surface of a semi-sphere 
in Fourier space. When we combine the exponential functions that contain the components of ro, we can identify the 
two-dimensional Dirac delta function 


J[dx D dy D exp (ix D (k m p - k Dx ) + iy D (k m q - k Dy )) = ( 2 vr ) 2 d(k m p - fc Dx , k m q - k Dy )). 


(5.14) 
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We can now evaluate the integral over p and q. 


maokn 


^B,0o( k c) = J^y /2 JJ dpdq M F ( km ( S ~ S °^ ' ex P ■ 5(k m p - fe Dx , k m q - fe Dy )) 


Ub,4>o ( k D) = 


t^B ,</> 0 ( k D) = 


iirao 


(2Tr) 1 / 2 k n 


dpdq T7FZ IT F (km [s (p, q) - So]) • exp (ik m M(p, q)l D ) -dip- q - 


z7rao 


exp 


M(p, q) 

(ib ^m-^Dx-^Dy) 


/CDx 

lT ,q ~~k' 

n, m A/ m 


(5.15) 

(5.16) 


(27T) 1 / 2 


1.2 _ 1,2 _ b .2 

K ra "'Dx ^Dy 


F ( fe Dx - fcmP0,fc D y - k m q 0 ,J A: 2 , - - fe£ y - fcmy 1 - Pq “ <?0 


(5.17) 


We used the identity of the delta function S(Ap ) = p^h(p), and inserted the definition of M = \/l — p 2 — q 2 and 
s = ( p,q,M). 


Ol 


to 

00 
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We may write the short form of equation 5.17 by setting 

&Dx = km'P and kv y = k m q 


fr n \ exp(ik m Ml D ) 

( k D ) = 7TiTlTW- Z 17- F ( k D - k m So) 


(2tt)V 2 k m M 
When we solve this equation for F , we get [3, 4 ] 11 


__ 12 7 Mk _ _ 

F (k D - /c m s 0 ) = - \ - -—exp ^ (k D ). 

V 7T ao 


complex factor 


(5.18) 


(5.19) 


5.2. Interpretation of the Fourier Diffraction Theorem 

Equation 5.19 shows a direct relation between the Fourier transform of the measured wave (Born 
approximation) C/b^o^d) and the Fourier transform of the inhomogeneity of the sample F(k). 

5.2.1. Surface of a Semi-Sphere in Fourier Space 

A closer look at the arguments reveals that the information defined by C/b^^d) is distributed 
along a semi-spherical surface with radius k m in Fourier space that is shifted by — fc m So, as is 
explained in the following. 

The spherical surface is defined by the wave vector with the magnitude k m . 

r .2 — n 2 _i_ h ,2 _j_ ;,2 

~ ft Dx < ft Dy T «-Dz 

Because k Bx and k By are given by the size of the detector image, we are left with a restriction for 
&Dz- This restriction forces the data FB,(/>o( k D) to be placed on a spherical surface with radius k m 
in Fourier space. Thus, the two-dimensional Fourier transform is projected onto a semi-spherical 
surface as depicted in figure 5.1. 

The shift in Fourier space is defined by the argument of F (ko — fc m So). The information of the 
Fourier transform of WB( r ) at the detector C/b^d) is shifted in Fourier space in the direction —so- 
The vector so is constant for a fixed value of <j >o and describes the propagation direction of the 
incident wave onto the sample. The magnitude of the shift is equal to k m which is the radius of 
the spherical surface described above. 


11 If the non-unitary Fourier transform and the sign of the Green’s function is taken into account, the term computed 
by Devaney et al. [3] (equation 48) differs by a factor of l/(k^ a ). This factor k ^ comes from a different definition 
of the data: /(r) = —k^ a O( r) in that paper. 
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Figure 5.1, Fourier diffraction theorem. The data C/b, 0 o (kyj ) (green) are projected onto a semi-sphere 
in Fourier space according to k ^ + k^ y + /cp z . The radius of the sphere is k m . The surface of the 

sphere is oriented along the s 0 axis in the case of coaxial illumination (detector aligned with incoming wave 
rto( r ))- The axes of the two-dimensional Fourier transformed detector image are labeled fcox and kn y . If the 
sample is rotated about one axis (red, here it is the y-axis) and the resolution of the setup is high enough 
(see sampling considerations below), then one obtains a horn torus-like shape in Fourier space (fig. 5.2). 


5.2.2. Sampling Considerations 

Diffraction tomography inherently increases and limits the resolution at the same time. As a 
result of the data distribution on semi-spherical surfaces in Fourier space, the maximum value 
of | k| = |ko — fc m so| is V^m- Apparently, the resolution is increased by a factor of \/2 when 
compared to the data measured. However, at the same time a maximum possible resolution is 
enforced by the restriction |k| < y/2k m in Fourier space. For the 2D case, the difference in resolution 
between projection tomography and diffraction tomography (maximum frequency in Fourier space) 
is visualized in table 4.1. 

It follows that all data distributed in Fourier space are located within a sphere of radius \/2k m 
according to 


^Dx + ^Dy < &m- (5.20) 

This frequency-limit in Fourier space infers a resolution limit in real space. The maximum frequency 
in Fourier space computes to 

f — /b ^’ m — (r 01 \ 

/max — v 2 ■ ^ — . (5.21) 

In other words, the optical resolution is limited to and depends on the refractive index of 

the surrounding medium n m . 

5.2.3. Comparison to two-dimensional diffraction tomography 

When comparing equations 4.24 and 5.19, it turns out that the Fourier diffraction theorem in two 
dimensions is similar to the Fourier diffraction theorem in three dimensions. Table 5.1 compares 
these two equations with respect to their dimensionality. 
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n k) = - 

2D (equation 4.24) 

4 ' \/^“^ B ^o( k D) 

3D (equation 5.19) 

Sinogram 

Fourier transform of complex 
scattered waves t/B,bo(^’Dx) 

Fourier transform of complex 
scattered waves f7B,bo(^Dx’ &D y ) 

Factor 

^ 2f& 

A —_ 

Iri M 

pxrff ih MItA 

A 

/i — 

C 

to 


M = W\A™ - fe Dx 

Aim 

M = TT J fc m - 4 - k Dy 

Ah m v 

Coordinates 

k = (/C X , ^z) 

= &Dx 

k = (/c x , fcy, /c z ) 

/p x — ^Dx^ — ^Dy 

k at ())q = 0 

h = \Jk 2, - k^ x - k m 

(semicircular arc) 

b, — lh.2 _ h.2 _ h.2 _ h, 

Aiz — y Ai m ^Dx ^Dy 

(semispherical surface) 


Table 5.1, Comparison of 2D and 3D diffraction tomography. Above the table is the generalized 
Fourier diffraction theorem for 2D and 3D. The object /(r) is rotated about the y-axis, which is the axis 
pointing away from the 2D ( x , z)-plane. The shape of the diffraction tomography formula is identical in 2D 
and in 3D. The only differences come from the different numbers of dimensions. For a comparison to the 
Fourier slice theorem, see table 4.1. 


5.2.4. Artifacts from Uniaxial Rotation 

Rotating the sample only about one axis results in missing-angle artifacts for 3D diffraction to¬ 
mography. According to the Fourier diffraction theorem, the resulting reciprocal volume in Fourier 
space has a shape of a missing apple core as depicted in figure 5.2 [23]. This results in directional 
blurring along the axis of rotation (y). In order to address directional blurring, regularization 
techniques need to be applied, as is described in e.g. [24-27]. 

5.3. The Limes of the Rytov Approximation 

In this section we show that diffraction tomography with the Rytov approximation converges to 
classical tomography when the wavelength A becomes small [2]. We start with the Fourier diffraction 
theorem which can be applied in combination with the Rytov approximation by calculating iiB(r) 
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Figure 5.2, Horn torus. As a result of the Fourier diffraction theorem, images acquired perpendicularly 
to the sample-rotating axis ( y ) fill a horn torus-like shape in Fourier space. The resulting reconstruction 
exhibits missing angle artifacts, as shown in figures 5.4 and 5.5 


from the measured scattered field it s (r) ~ %(r) using equation 3.39. 


ub{ r) = rto(r)ln 


ao f ttR,(r) 
V a L M o(r) 

= M 0 (r)<^R(r) 


+ 1 


The Fourier transform of the equation above at the detector r — > ro yields 

1 


^B,0o(k d) = JJd z r D u o (r)^R,0 o (r D ) exp(-zk D r D ) 

= J 0 (rD)exp(-i(k D - k m s 0 )r D ) 
= «o &R,cpo (^d — k m s 0 ) 


(3.39) 

(3.40) 

(5.22) 

(5.23) 

(5.24) 


where (jr^^d) is the two-dimensional Fourier transform of ^R, 0 o (rD)- For the Rytov approxi¬ 
mation, equation 5.19 then reads 12 


F (k D - k m s 0 ) 



iMk m exp(-ik m Ml D )p R ^ 0 (k D 


fcm®o) • 


(5.25) 


By replacing k^j = kn — & m so one can rewrite this equation such that the argument of F( k) does 
not anymore contain the shift in Fourier space. 


F (k' D ) = - exp(-ifc m M*/ D )^R j( £ 0 (k' D ) (5.26) 

M* = ^l-tp' + PoY-W + qv) 2 (5.27) 

k' D = k m (p',q',M') (5.28) 

= k m (p- p 0 ,q- qo,M - M 0 ) (5.29) 


12 Note that kn is in a plane as ro is defined in the detector plane. Therefore, (kn — fc m so) is a two-dimensional 
Fourier transform. However, F (ko — k m so) is three-dimensional because it contains the fc z component. 
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The data are still distributed on the surface of a semi-sphere in Fourier space. This information is 
hidden in the unintuitive restraint from M* 

C p'+Po) 2 + (q' + qo) 2 <l • (5.30) 

This expression can be simplified by noting that for = 0, is parallel to so - our measuring 
^-direction (Mo = 1) and therefore Po = Qo = 0 and 

p + q <1 

M* = y/l - p ,2 - q’ 2 
ko = k m (p, q, M - 1) 

Note that the magnitude of the vector k'j-, is not fe m 

|k' D | = 2fc m (l - M) + fe m (5.34) 

and thus the restrictions for the k z component of k^ are different. It is shifted by —k m in the fez- 
direction. The shift in Fourier space is not obvious, but it is included in the restrictions for the fe z 
component of k'j-,. 

In the short wavelength limit A approaches zero and thus fe m , the radius of the semi-spherical 
surfaces in Fourier space, goes to infinity. The result is a data distribution that is identical to that 
of the Fourier slice theorem. Furthermore, short wavelengths imply fe x , fe y <C fe m and we may write 
M* —> 1 [2] to obtain 


F(kb) 

= “ \f^ ikm ex P(-^mZD)^R,^ 0 ( k D) 

(5.35) 

^R,*j( k D) 

= 2tt ff dxD(lyD o( r D)exp(-7;k / D r D ) 

(5.36) 

rfki,) 

= ( 2vr ) 3/2 /[ j dxd y (iz /( r ) exp( *k D r). 

(5.37) 


If we perform a two-dimensional Fourier transform of F (k^-,) in the reciprocal detector plane after 
inserting 5.36 into 5.35, we get a projection along the z-axis smeared out by a periodic exponential. 

ik m exp(-?:fe m / D )^ R ^ 0 (r D ) (5.38) 

Where the interval [—A, A] is the domain of /(r) along the ^-direction. Here, we used fe’Dz = 
M — 1 —> 0. The exponential factor on the right hand side only yields a phase offset in the detector 
plane. We may set Zd to zero and arrive at [2] 

r-\-A 

/ dzf(r) = — 2ife m ( / ? Ri 0 o (r D ). (5.39) 

■J-A 




(5.31) 

(5.32) 

(5.33) 
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The right hand side contains the complex Rytov phase in the detector plane v ? R,0 o ( i ‘d) = ?< ^( r D) 
(see section 3.2). The left hand side contains the scattering potential from equation 2.9. 


/(r) = kl 


= K 


n(r) 


n i, 


- 1 


1 + 


£»r 


Wn 


- 1 


en(r) 2 <€„(r) 2 k"L 


n r 


e n( r ) 


(2.9) 

(5.40) 

(5.41) 


Where e n (r) is the local refractive index variation from the surrounding medium n m . By writing 
the right hand side of equation 5.39 as an integral over d<h(rD), we thus get 

/ +A dz r 

2kl 1 e n (r) — = / 2A: m d$(r D ) (5.42) 

-A n m J 

d z 

2^m e n( r ) — = 2& m d<f>(r D ) (5.43) 

€ 4±dz=U*(r»). (5.44) 

A 27T 

This relation describes the phase change that occurs over a distance dz, as derived in section 3.1, 
equation 3.12. Thus, in the limit of high k m , the Fourier diffraction theorem with the Rytov 
approximation becomes the Fourier slice theorem, which requires that the data (fd$( ro)) recorded 
at the detector ro are computed from line integrals ( fdz ) through the sample e n (r). 


5.4. Derivation with the Fourier Transform Approach 

An alternative way to derive the Fourier diffraction theorem (equation 5.19) is the convolution 
approach in Fourier space. This approach uses the convolution theorem for Fourier transforms. We 
only briefly discuss this approach. For a thorough discussion, see [22]. 

tt B (r) = Jid 3 r' G(r - r')f(r')u 0 (r') = (G * f u 0 )( r) (5.45) 

Equation 5.45 describes the Born approximation '«b(i") as a convolution (*) of G(r) with f(r)uo(r). 
In Fourier space, we may write Ub (k) as 13 

C/ B (k) = (2vr) 3 / 2 G(k) • CfS)(k), (5.46) 

where (/ito)(k) denotes the Fourier transform of /(r)ito(r), and G(k) is the Fourier transform of 
G( r). We find that 

(fu 0 )(k) = ^ 9 ^° 3/2 Jd 3 r'f(r') exp (-i [k - k m s 0 ] r) = a 0 F( k - k m s 0 ) (5-47) 

13 The factor (27r) 3 ' 2 originates from the unitary angular frequency Fourier transform. The non-unitary angular 
frequency and the unitary ordinary frequency Fourier transforms do not show this factor. 
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is simply a shifted version of the Fourier transform of the inhomogeneity /(r). In order to find 
G(k), we calculate the Fourier transform of the inhomogeneous Helmholtz equation for the Green’s 
function (equation 2.10). 


^G(k) 


( 2 7 t) 3 /2 /^ r ' V2G ( r/ ) eX P( ?:kr ') + (^ 3/2 j d3r ' k m G ( r ') exp(ikr') = 


( 2 7t) 3 / 2 


d 3 r' 5 (r) exp(ikr') 


(5.48) 


Using integration by parts for each Cartesian coordinate and by considering the asymptotic behavior 

exp (ikr) exp(ifc|x|) 


as x 


oo 


47T7’ 47T \ x \ 

one can show that the first integral in the equation above becomes [22] 

1 Ob' V7- 


d 3 r' V^G(r') exp^kr 7 ) = — ArG(k). 


( 2 7t) 3 / 2 

Thus, we obtain the Fourier transform of the Green’s function 

G( k) = 1 1 1 


1 


( 2 7t) 3 / 2 k 2 - k ^ ( 2 7t) 3 / 2 k 2 + k 2 + k 2 - k^ 

We combine the equations above and obtain 

fr B (k)= Q ° -Ok-fcmSo). 

/€ n yv. 


(5.49) 


(5.50) 


(5.51) 


(5.52) 


This equation looks very much like equation 5.19. However, here the entire held UB(k) is on the 
left hand side of the equation. Since we measure the held at the detector t/B )( /, 0 (kD), we perform 
an inverse Fourier transform to normal space, then take the held ub (r) at the detector, and finally 
back-transform to two-dimensional Fourier space. We rotate our coordinate system such that so 
matches the £i z -axis. The integral over dk z can be evaluated using the residue theorem. We integrate 
around the singularity at k z = yj, A; 2 , — k 2 — k 2 . 


u B (r) = 


«o 


( 2 7t) 3 / 2 

2 7riao 

( 2 7t) 3 / 2 


dk x dkydk z ^ ^ j ^ ^m®o) 


k 2 - kl 


dk x dk y - 


expl i 


Is* T* I 'll I / A <y 

ruyy \ 


9 / Z-2 — k 2 — k 2 

Zj \ / rvy 


F (k - A; m s 0 ) 


(5.53) 

(5.54) 
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Here, we used the following partial fraction expansion to obtain the singularities. 

1 _ 1 / 1 1 

i,2 _ y2 ~ ou l u _ y i, y 

''in z/ ’in \">z ''m l^z i ''in 

k'm = — ky (5.56) 

We identify k as ko = (.td , 2/d , zd = Zd) (at (/>o = 0) and perform the two-dimensional Fourier 
transform at the detector. We obtain delta functions for x and y that we use to solve the integral 
for k x and k y 14 . With M = y— k^ x — k^ y we arrive at equation 5.19. 

f T n x 0 exp (ik m Ml D )~ 

C/B^ 0 (k D ) = 2vr (27r)1/2 _ 2 -^-F(k D - fe m s 0 ) (5.57) 


(5.55) 


5.5. Backpropagation Algorithm 


In this section we derive the three-dimensional analog to the two-dimensional backpropagation al¬ 
gorithm as described by Devaney [3] and section 4.2 in this script. Our goal is to solve equation 5.19 
for /(r) such that we do not need any interpolation in Fourier space. In order to do that, we need 
to perform a change of coordinates. The coordinate system has its origin on the rotating axis y. 
We define our new coordinate system in the Fourier domain analogous to the two-dimensional case 
in three steps. 

First, we construct our coordinate system such that we can eliminate any reciprocal distances and 
express any position of the recorded data in terms of angles. This is possible, because the Fourier 
diffraction theorem states that every point in Fourier space is placed on a semi-sphere with radius 
k m centered at — fc m so. Second, we define the projection angle (fro as the angle that corresponds to 
the angular sample position from which the projection was recorded. The angle <fro coincides with 
the direction of the incoming plane wave So- Note that (fro is a two-dimensional angle. However, 
if the sample is rotated only about one axis, the angle (fro will become one-dimensional. We only 
consider the rotation about the y-axis as shown in figure 5.3a. Figure 5.3b depicts the two angles 
which define the semi-spherical surface in spherical coordinates relative to the rotated coordinate 
system ((fro)- Third, we set the polar angle as 9 and the azimuthal angle as ifr. For (fro = 0, the 
non-shifted semi-sphere lies within the half-space k z > 0. The polar angle 6 is then measured from 
the positive k z - axis and the azimuthal angle ifr is measured from the k x -ax is in positive direction of 
ky. The corresponding variable domains are 


9 G 


7T 7T1 



ifr G [0, tt] . 


Let us first write the rotation matrix D 0() for the rotation of the sample about the y-axis. 


D<t>o 


cos <fro 0 — sin (fro 

0 1 0 

sin <fro 0 cos (fro 


(5.58) 


14 


Factors: the two-dimensional Fourier transform adds a factor of l/(27i) and the delta functions contribute with 27r 
each. 
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in the rotated system 

Figure 5.3, Coordinate transforms for 3D backpropagation. a) The coordinate system of the detector 
(red) is rotated about the y -axis which corresponds to a rotation about the fc y -axis in Fourier space, b) The 
position on the semi-spherical surface with the radius k m is described in a spherical coordinate system with 
the polar angle 0 and the azimuthal angle ip. 


The matrix D ( . H] rotates a vector s' counter-clockwise through the angle (p o about the fc y -axis. The 
vector that describes the senri-spherical surface in spherical coordinates is given by 

s' = (sin0cos?/’, sin?/’sin 9, cosd) T (5.59) 

where () T denotes the transpose of the vector, i.e. s' is a column vector. We transform this vector 
into the rotated coordinate system by calculation of 

( sin 6 cos (/> o cos ?/> — sin </>o cos 6 \ 

sin ip sin 9 1 . (5.60) 

sin cpo sin 9 cos ip + cos <po cos 9 J 

We use the change of coordinates in Fourier space to rewrite the inverse Fourier integral for /(r) 
(eq. 5.19) 

^ / 2 x]\/[]^ / ,_ 

F{ k D - k m s 0 ) = - J -- exp(-ifc m Af'/ D )C/ B (kD) (5.19) 

V vr a 0 

/( r ) = ^)3/2 /// d * k F ex P( ikr ) ( 5 - 61 ) 

where we used 


k = k m ( s - s 0 ) 

(5.62) 

k D = k m s 

(5.63) 

s 0 = (-sin cpo, 0, cos(/> 0 ) T . 

(5.64) 
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The integral over d 3 k can be expressed as an integral over (p o, 9, and ip by means of the Jacobian 
matrix J. 


dk x dk y dk z = |det(J)| d<pod9dip 
dk x dk y dk z 
8 (pod 9 dil’ 

The Jacobian matrix is the matrix of all first-order partial derivatives of the vector 


( kx \ 


f sin 9 cos (po cos ip — sin (po cos 9 + sin (po \ 

k y 

— ^m(s Sq) — k m 

sin ip sin 9 j 

V kz ) 


^ sin (po sin 9 cos ip + cos (po cos 9 — cos (po / 


and computes to (column wise): 


dk x dk y 8k z f 

Wo = " 


dk Y dkydk z 

— rC n 


89 


dk x dk y 8k z 

— fCrI 


dip 


— sin <pQ sin 9 cos ip — cos (po cos 9 + cos (po 

0 

sin 9 cos (po cos ip — sin (po cos 9 + sin <po 

cos (po cos ip cos 9 + sin (po sin 6 
sin 'ip cos 9 

sin (po cos ip cos 9 — sin 9 cos (po 

— sin ip sin 9 cos (po 

sin 9 cos ip 

— sin (po sin ip sin 9 


(5.65) 

(5.66) 

(5.67) 

(5.68) 

(5.69) 

(5.70) 


For the integration by substitution we thus need to consider the following factor for the integral 15 

|det(J)| = |— /^(sin#) 2 cos -01 (5-71) 

and we may rewrite the integral for /(r) as 

1 1 f 2n Z* 71 " / i +7r/2 ^ 

/( r ) = 0 /0 \ 3/2 / d( ^ 0 d0 / I (sin 6») 2 cos V'I F(k)exp(ikr) (5.72) 

2 ( 2 ? t) 6 / 2 Jo Jo J-n/2 

where we have introduced a factor of 1/2 to compensate for the double-</>o coverage in Fourier 
space as we rotate the semi-sphere about the k y axis from 0 to 27r. Note that this factor only 
holds for equidistant coverage over 27T. In other cases, weighting factors need to be introduced for 
each angle (po■ We have to perform one more coordinate transform to arrive at an integral over 
the coordinates fco x and &d y - These are the coordinates in Fourier space that correspond to direct 
Fourier transforms of the two-dimensional recorded images at angles (po■ This coordinate transform 
is also depicted in figure 5.3 and reads 


kf) x = k m cos ip sin 9 (5.73) 

kny = km sin ip sin 9. (5-74) 

15 This determinant can easily be calculated using SAGE (http://www.sagemath.org/). 
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For 0o = 0, the &Dx-axis coincides with k x and the y -axis with k y . The Jacobian determinant 
Jd for this coordinate transform at the detector then computes to 


11 , / J \ I 3/C£) x 02C£)y 

|clet(jD)l= MW 

(5.75) 

=AinJ cos#sin# 

(5.76) 

06 [ “" / i’ + " / 2] ^ COS #| si n#|. 

(5.77) 


In the second variables change of integral for /(r) we identify k m cos8 = k m M. 

1 


dOdJj = 


rdk^dk^y 


fc^cos#! sin#| 
klM\ sine\ dkBxdkBy 


(5.78) 

(5.79) 


By applying these transforms, one obtains an integral for /(r) over the rotation angle 0o and the 
coordinates of the Fourier-transformed image fco x and /cd y - 


/(r) = J 1 


2 (2vr) 3 / 2 J 0 


2-7r rkm rk m |L I ^ 

c20 o / (2 /cdx / d/c Dy ^^-F(k)exp(2kr) 


'-fcn 


' -fcn 


M 


(5.80) 


The /^-component of the vector ko can be expressed by &Dx> ^D y) and 0o- We express k in terms 
of ko and sq, which can be separated into lateral (tj_) and axial (sq) components. 


k = 


t± = 


so = 

So • t_L = 


kD — 2c m so = k m (s- s 0 ) 
2cDxt± + km(M - 1) So 

^COS0O, sin 00^ 

(-sin0o, 0, cos0o) T 
0 


(5.81) 

(5.82) 

(5.83) 

(5.84) 

(5.85) 
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By inserting the relations above, we may may express /(r) as 


1 1 f 2rr f km f km Ifcnl — 

/( r) = - 3/2 / dt(t> o dk Dx dk Dy -Tf~ F (7 Dx ty + k m (M - 1) s 0 ) exp[i(£: Dx t ± + k m (M - 1) s 0 )r], (5.86) 

4o J fc m J fc m 

Now we insert equation 5.19. Note that the two-dimensional function t/B, 0 o (kD) was measured in the detector plane and 
thus we write Up,^ 0 {kp, x , kp> y ). Equation 5.19 describes the Fourier-space distribution of the measurement data from one 
single projection at an angle 4> o- 


^ / 2 iMk m . ^ 

1 (h i) fc m s 0 ) — W exp( ik m j\11 \))17b, Po (k \) x , /i’Tjy j 

V vr ao 


(5.19) 


The integral over <ji>o depends on the recorded data Up J;( p 0 . Furthermore, our choice of coordinates requires the exponential 
factor in the Fourier transform to be dependent on (f> q. 


/(r) = 


—ik n _ 

( 2 vr) 2 a 0 J 0 


r*27r 


# 0 / d/c Dx / d/c Dy |/c D x|%,(/>o(^Dx,^Dy)exp(-z/c m M/ D )exp[i(A;Dxt_L + A: m (M- l)s 0 )r] (5.87) 


'-Ain 


' Ain 


with 

k m M = \Jk^ — k^ x — (5.88) u, 

Equation 5.87 is the three-dimensional backpropagation algorithm of uniaxially rotated samples. An application of the 
3D backpropagation algorithm is illustrated in figure 5.4 for a centered sphere from computations that are based on Mie 
theory. Note that, as could be seen for the 2D case in figure 4.1, the reconstruction with the Born approximation breaks 
down and that the reconstruction quality with the Rytov approximation is dependent on the total number of projections 
that were acquired for the sinogram. Details of the simulation and reconstruction process can be found in figure 5.5. 
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5 Three-Dimensional Diffraction Tomography 






(d) Born x = 0 

(200 projections) 


(e) Rytov x = 0 
(50 projections) 


(f) Rytov * = 0 
(20 projections) 





(g) Born y = 0 

(200 projections) 


(h) Rytov y = 0 
(50 projections) 


(i) Rytov y = 0 

(200 projections) 


Figure 5.4, 3D Backpropagation. Image (a) shows the refractive index distribution of a dielectric sphere. 
White indicates the refractive index of the medium and black the higher refractive index of the sphere. The 
forward-scattering process was computed according to Mie theory. Amplitude and phase of the scattered 
near-field are shown in (b) and (c). Because the sphere is rotated about its center, the amplitude and phase 
are identical for every slice of the sinogram. Images (d), (e), and (f) show the reconstruction of the refractive 
index distribution using the Born and Rytov approximations in the x = 0 cross-section. Images (g), (h), 
and (i) show the corresponding reconstructions at the y = 0 cross-section. The reconstruction in (e) and (h) 
were performed with 50 projections only. See figure 5.5 for details. 
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5 Three-Dimensional Diffraction Tomography 
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Figure 5.5, 3D Backpropagation cross-sections. The refractive index of the medium is n m = 1.0 and 
the local variation inside the sphere is e n (r) = n(r) — n m = 0.006. The radius of the sphere is 14A (vacuum 
wavelength A). The scattered wave is computed at an optical distance of Zd = 20A from the center of the 
cylinder and sampled at A/3.1125 over 250 pixels using the software GMM-FIELD [28]. The refractive index 
map is reconstructed on a grid of 250 x 250 pixels from 200 projections over 27r using the software ODTbrain 
[21]. The reconstruction with the Born approximation (not shown) follows the behavior as seen in figure 4.2. 
The 2D cross-sections are shown in figure 5.4. Artifacts due to the incomplete 27r-coverage are visible along 
the y- axis, as briefly discussed in section 5.2.4. 
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6 Implementation 


5.5.1. Comparison to the Two-dimensional Backpropagation Algorithm 

Table 5.2 summarizes the similarities and differences between the 2D and the 3D backpropagation 
algorithms with respect to the dimension of the problem. 



' w ~£,G H 

2D (equation 4.41) 

r 2n ^ 

/ 00o exp (iB) |fc Dx | (ko) 

Jo 

3D (equation 5.87) 

Sinogram 

^o^d) 

ID Fourier transform 
of complex scattered waves 

^B,0 o (^Dx) 

2D Fourier transform 
of complex scattered waves 

&Dy) 

Integral 
dK ' D 

(/ dA °) “ V2^I dkDx 

( I ^ D ) = 2tt f J ^ Dx ^ Dy 

Exponent 

B 

B — — k m AflY) + 

M = TT\/ k l~ k D x 

'*T11 

/c Dx t±r + k m (M - l)s 0 r 

M = ~T~ \J k m “ ^Dx “ ^Dy 

a, m v 

Vectors 
r, s 0 , t ± 

r = (x, z) 

s 0 = (— sin0o, cos 0o) 
t_L = (cos 0o, sin 0 O ) 

r = (x,y,z) 

s 0 = (-sin0o, 0, cos0 o ) 
t_L = ^COS0O, sin00^ 


Table 5.2, Comparison of 2D and 3D backpropagation. The table compares the 2D and 3D back- 
propagation algorithms. The generalized backpropagation formula above the table shows that the structure 
of the algorithms are similar in 2D and 3D, as previously noted in table 5.1. The only differences originate 
from the different numbers of dimensions. See table 4.2 for a comparison between the backpropagation 
algorithm and the backprojection algorithm. 


6. Implementation 

The actual algorithms for backprojection and backpropagation use a trick for the reconstruc¬ 
tion. We do not need to numerically integrate the equations for the backprojection (1.16), the 
two-dimensional (4.41), and the three-dimensional (5.87) backpropagation algorithms. Instead, 
we identify the Fourier transform for the reciprocal vector ko and perform the reconstruction 
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6 Implementation 

projection-wise for each rotational position <p q. In this section, we subsequently derive the actual 
reconstruction algorithms from the equations of the previous sections. 


6.1. Backprojection 

The backprojection algorithm, as introduced in section 1, inverts the Radon transform, i.e. it is 
used to reconstruct an object from equidistant projections. Figure 1.2 depicts the process from 
image acquisition to image reconstruction with the backprojection algorithm. The backprojection 
algorithm is based on the Fourier slice theorem [7, 9] and runs as follows: 

1 . Each projection is filtered with a ramp filter (|A:dx|)- 

2. The filtered projection is backprojected onto the image volume according to its acquisition 
angle. 

3. The sum of all backprojections constitutes the reconstructed image. 


We start from the formula for the inverse Radon transform. 


f(x,z) 


1 

27T 



d<t >0 |&Dx 


P<j >o (^Dx) 

\/27T 


exp[ifeDx(a: cos fio + z sin c^o)] 


(1.16) 


The data in real space at r = (x,y) are computed from integrals over kr> x and <f> o- We can 
introduce a coordinate transform D_ ( j >0 that rotates r through the angle —fi o along the y-axis, 
such that x<f, 0 = x cos fio + z sin fio . In the resulting equation we identify a one-dimensional inverse 
Fourier transform 


f(x,z ) = 


2tt 


#o Jdk Dx |/cdx| exp[ifc Dx a^ 0 ] j • 


FFT^{|fc Dx |P 0o (fc Dx )} 


( 6 . 1 ) 


We have effectively replaced the integral over k-£> x by a one-dimensional inverse fast Fourier trans¬ 
form (FFT^) and a rotation in real space (H_ 0 O ). We replace the remaining integral over (f> o by 
a discrete sum over N equidistant projections and obtain 


f(x,z) 


1 


N 

^><(>0 D. 
3-1 


^{fft^ { |/c D x | (fc Dx )}} 


( 6 . 2 ) 


with the discrete angular distance A^o = vr/A and the discrete angles <f>j = j- A0o (j = 1,2,..., N). 
A numerical method that implements equation 6.2 is much faster than the direct computation 
of equation 1.16, because it can make use of the fast Fourier transform. The common name 
filtered backprojection algorithm comes from an interpretation of equation 6.2. The one-dimensional 
detector data are ramp-filtered and then projected over the entire two-dimensional reconstruction 
plane according to the angle they were acquired. 
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A Note on the Ramp Filter 

The ramp filter |&Dx|, derived in section 1.3 can also be motivated by the point spread function of 
the backprojection process, which is t/|r|. By filtering each projection with the ramp filter |/cdx|> the 
effect of the point spread function can be reversed. To illustrate this, let us consider a reconstructed 
image that has not been filtered (fo). Then we may write a formula for the filtered image / by 
means of the convolution 

/o = /*n- (6-3) 

l r l 

Inverting this formula using the Fourier transform J 7 , leads to 


fo 


f 



T~ x CF(/o) • |k| /2tt) 



(6.4) 


(6.5) 


The derivation shows, that the image in real space / can be computed from the non-filtered back- 
projection /o by multiplication with a two-dimensional ramp filter |k|. Filtering each projection 
with |A;dx| is mathematically equivalent and yields better results when working with discrete data 
sets. With regard to detector symmetries or sample type, other filters (e.g. cosine, Hamming, 
Hann) have been developed. These filters are usually heuristic and have different effects on image 
contrast or noise. 


6.2. Two-dimensional Backpropagation 

The backpropagation algorithm in two dimensions (r = (x,z)) can be derived analogous to the 
backprojection algorithm. An additional difficulty is introduced by the fact that the projection is 
now multiplied by a phase factor which is dependent on the distance from the center of the object. 
Let us consider equation 4.41. 


#o |&Dx| U Bj(j)0 (kn x )exp(-ik m Mln)x 
x exp[i(& Dx t_L + k m (M - 1) s 0 )r] 


(4.41) 


t_L = (cos 0o, sin 0o) 
s 0 = (-sin0o, cos0 o ) 
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6 Implementation 

We begin by introducing the rotation -D_0 O through —cf >o along the y-axis that transforms r to r^ 0 . 


T <h 

( X <j>0 ’ Z <t>0 ) 

(6.6) 

X <f>Q 

= x cos <f>o + z sin 

(6.7) 

Z <!>0 

= —x sin <f>o + z cos rf>o 

(6.8) 

t_L • r 

x 4>o 

(6.9) 

So • r 

= Ztj, o 

(6.10) 


/( r ) = 


ik m 

ao(27r) 3 / 2 




'dk Dx |&Dx| (&Dx) exp(-ifc m M/ D ) exp[i(fc D x®</> 0 + &m(M - 1 )z <j>0 )\ 


( 6 . 11 ) 


Because of the factor k m (M — 1) in the integral, we cannot proceed exactly as we did in the 
previous section. The inverse Fourier transform needs to be performed for all coordinates before 
the rotation. 


/( r ) 


27t • a o 



x £>_ 


0 o {fFT 1 p ||fc Dx | C/ B ,</> 0 (^Dx) exp(—exp[ifc m (M - 1)^ 0 ]}} 


( 6 . 12 ) 


In practice, the measured scattered field u s (r), which is approximated by the Born approximation 
«B(r), is already background corrected. We insert the incident plane wave uo(/d) = a o exp(ifc m /n) 
and thus, obtain an additional factor exp(+ik m li)) in the Fourier filter. 


^(fex) FFTid{ub(®)} 


«0 


ao 


= FFT 1D 


mb(T) 
u o(Id) 


■ exp (ik m ln) > = 


^B(fcDx) • exp(zfc m ?p) 
uq(Id) 


Using this substitution and the discretization of the integral over according to the previous 
section, we obtain the backpropagation algorithm in two dimensions 


N 


/«—„D 

3 =1 


FFT 


ID (I^Dxl exp [ik m (M - 1) • (Zfa - Z D )] 


(6.13) 


with the discrete angular distance A0o = 2n/N and the discrete angles 4>j = j ■ Afo (j = 1,2,..., N ) 
When comparing this equation with the backprojection algorithm from equation 6.2, we can see 
one major difference besides the different filter. The inverse Fourier transform needs to be calcu¬ 
lated separately for every distance z$.. In practice, one first calculates the one-dimensional signal 
UB^jik^x)/ uo(Ib) ■ exp [—ik m (M — 1)/d] and then expands the signal by one dimension through 
the application of the second filter exp[ifc m (Af — 1 )^.]. The inverse Fourier transform is then 
computed along the axis with constant z ( p j . The resulting two-dimensional data are rotated by (j)j 
and added to the reconstruction plane. The name backpropagation comes from an interpretation of 
the z ( j )j -exponential, which looks like a propagation in z^ -direction. 
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6.3. Three-dimensional Backpropagation 


With the previously described two-dimensional algorithm, it is now straight-forward to derive the 
three-dimensional analog. We start again from the integral representation. 


/( r ) = 


— ikr. 


p 27r rh m rk m 

d(j ),o / d/cDx / dk Dy |fc Dx | x 

•' k m k m 


(2tt) 2 a 0 Jo 

X Hb,</, 0 (£;dx, ^’Dy) exp(-i/c m M/ D ) exp[i(/c Dx t± + fc m (M - 1) s 0 )r] 


(5.87) 


t ± = (cos^o, sin (/) 0 

V KDx 

s 0 = (— sin fa, 0, cos ^o) 


T 


T 


In this section, we only consider the rotation of the sample along the y-axis. This introduces 
artifacts, as briefly discussed in section 5.2.4. The rotation through — </>o is described by the rotation 
operator D_^ 0 that transforms r^ 0 to r. 


r 

= (ac, 2/, ^) T 

(6.14) 

r 0o 

(T0o > 2/0o) z 4>o ) 

(6.15) 

x 4>o 

= x cos (f>o + z sin 4>o 

(6.16) 


= y 

(6.17) 

^00 

= —x sin <f>o + z cos 4>q 

(6.18) 

t ± • r 

_ k]0y 

— x <po T t ' 2/0o 

«Dx 

(6.19) 

So • r 

= Zfj, o 

(6.20) 


r*27T 


/( r ) = 7^ / d( ^° D ~^ 
(2vr)-a 0 J 0 


ffcm rk m 

d^Dx / dfcoy |A;dx| x 


'-fcn 




xt7B,0o(^’Dx, ^Dy) exp(-i/c m Ai7 D ) exp[i(feDx®</. 0 + ^D y 2/0 O + /c m (M-l)^)]} (6.21) 

Here, we identify the two-dimensional Fourier transform over /cdx and /coy and simplify the equation. 


/(r) = 


ik n 

2vr • a 0 Jo 


r*27T 


dcj>o x 


x T>_0 o |fFT 2 q |Hb,</) 0 (A:dx) fe Dy ) |^dx| exp (~ik m Ml D ) exp[ifc m (M - 1)^ 0 ]|| (6.22) 
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Following the derivation in two dimensions, we identify the incident plane wave uo(/d) and discretize 
the ^-integral 

3 =1 

x D_J FFT~' | ^ B ’ 0 ^’ fcDy) |fc Dx | exp[ik m (M - 1) • (z^ - l D )] j j (6.23) 

with the discrete angular distance A0o = 2n/N and the discrete angles 4>j = j • A</>o (j = 1, N ). 
Equation 6.23 and equation 6.13 only differ in the number of dimensions. In the three-dimensional 
case, the acquired data are two-dimensional projections through the sample. The fast Fourier 
transform is used again to speed up the computation. The subsequent rotations are performed for 
three-dimensional volumes which are then added to the reconstruction volume in the sum over (pQ. 


Conclusions 

This manuscript reviewed the theory of diffraction tomography. We discussed the wave equation 
with the Born and R.ytov approximations. Next, we compared the mathematical formalism of two- 
dimensional diffraction tomography with the formalism of classical tomography and pointed out 
the major differences between them. Subsequently, we derived the full three-dimensional equivalent 
of diffraction tomography and discussed the implications for biological samples. Finally, we gave 
details on how to implement the backpropagation algorithm in software. This script is intended 
as a literature review, establishing the mathematical foundations of diffraction tomography with 
a coherent notation. The reconstruction methods described are applicable not only to optical 
diffraction tomography but also to ultrasonic diffraction tomography whose underlying principle is 
the scalar wave equation. In general, diffraction tomography opens new ways for three-dimensional 
tissue imaging and with growing computational power at hand, will become a feasible and flexible 
tool in modern imaging. 
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A. Notations in Literature 

A.l. Two-dimensional Diffraction Tomography 

Devaney and Slaney derived the Fourier diffraction theorem with a slightly different notation [3, 
18]. This section illustrates the main differences in notation and translates between them. Besides 
different variable names and a differing definition of the object function/scattering potential /(r), 
the Devaney and Slaney use the non-unitary angular frequency Fourier transform, whereas this 
script uses the unitary angular frequency Fourier transform. The Fourier diffraction theorem from 
equation 4.24 reads 


^ 2 ik ^ 

F(k m ( s - s 0 )) = - \ (fc Dx ) exp(—i/c m M/ D ). (4.24) 

V vr ao 

Using the notation of Slaney, we start with the Fourier transform of the scattered wave A ([18], 
eq. (62) 16 ). 

J A(x,y = l)e~ iax dx = f(a,l ) 

f (a',l) = J e~ ta ' x j F, (a, l)e iax dadx 

Note that Slaney uses the non-unitary angular frequency Fourier transform. With the identity of 
the delta function 

S(a - a') = 1 [ e’( Q - Q >cte 

2vr J 

we may write 


r (a',l) 


27 n 5 (a — a')Ti(o/, l)da = 27 rFi(a / , l). 


Using equation 60 from [18], we get 1 ' 

A (x,y = l)e~ lc 


I 


c dx = 


l fQ(«, - fc 0 ) 


2 V / ^r r 




([ 18 ]) 


Note that in Devaney et al. [3], equation (20) is equivalent to the equation above, when replacing 
O with —k 2 0 and performing a one-dimensional Fourier transform. Now we perform the following 


16 The factor 2 ' is missing in [18]: The derivation with the convolution theorem does not seem to be correct. Equation 
(47) in [18] should read V's(A) = G( A) * -77(A), where 77(A) is the Fourier transform of H(f) = {O(r)ipo(f)} . 
The additional factor of 2-7T seems to originate from a delta function in equation (50). The full three-dimensional 
Fourier transform approach is derived in section 5.4. 

1 'The exclamation mark (!) in this formula points at the sign of equation (60) and (61) in [18], which should be 
switched. This follows from the path integration described in the paper. 
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translations to get from Slaney’s notation to the notation in this script. 


detector position: 

(x,y = l) -a (. t d , zd = b) 
l ->• Id 

correlating factor: 

y/k% — a 2 —> k m M 
Fourier transformed object: 

0(a, y/k% — a 2 - k 0 ) -A 27 tF (k m ( s - s 0 )) 

(object function defined at [18], eq. (10)) 

[ ip s (x,y = l)e~ iax dx -> *—Ub^o^Dx) 

7 oo 


By replacing these things, we get the equivalent to equation 4.24. 


(27T) 1 / 2 

ao 


UB,<fr 0 (kDx) 


2mF(k m (s - s 0 )) 
2k m M 


exi)(ik m Mly)) 


The backpropagation formula (equation 4.41) can be translated via the exact same way. Using the 
previously mentioned publication by Devaney [3] and combining the equations (35), (36) and (29a) 
therein, the backpropagation formula reads 


0 L D(r) 


-i- f#, o ["(IK | K |r^ 0 ( K ,o;)e^+^- fe )^). 
J — 7T J — k ^-N/' 1 

/ oo 

-OO -V-' 


4 

n p iklo . . 

-Ufa= l o,u) 


kU 0 (uj) 


By inserting all of these translations we find the Fourier transform ujj^ (k. uj) of the scattered wave 

= lo,w). 


Old ( r ) = 


(2vr) 


d(j )o / dn |k| 


l-k 


kU 0 (u) 


d? y = l 0 ,u)e- iK ^ 




The Theory of Diffraction Tomography 


50 



A Notations in Literature 


We now perform the following translations. 


variables: 

lo 

k 

7 

factors (notation): 

OldW 


Id 

km 

k m M 

A/ m 


'4>o 

amplitude: 

U 0 (u) 

vectors: 

£ 

v 


(sign and factor k^ n defined at [3], eq. (12)) 
-A V2TrU B ^ 0 (k Dx ) 

(non-unitary to unitary Fourier transform) 


a o 

t± 

so 


We arrive at the previously derived backpropagation formula. 


/( r) = ~ 


ikr> 


r2n 


, n ,o /9 dk Dx / #o |fcDx| CW*Dx) exp(-ifc m MZ D )x 

Clo^TT) 6 / 2 J J o 

X exp[i(/c Dx t_L + k m (M - 1) s 0 )r] 


(4.41) 
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A.2. Three-dimensional Diffraction Tomography 

In section A.l we looked at two notations in literature that describe the two-dimensional case. The 
three-dimensional case is briefly discussed by e.g. Devaney [3]. We show the equivalence of the 
notations by performing a Fourier transform of equation 5.18. 


UB,<t> 0 ( r D) = 


2tt 


exp(?'k D r D ) 


z7rao 

(27r)V 2 k m M 


F (kn — k m s 0 ) 


The following translations have to be performed to translate between the notation of Devaney and 
ours. 


detector position: 

(z D , 2 /d, zd = Id) -a (x, z,y = Q 
Id A y l Q 

UB,<t > 0 (rD) -A t/g s) {x, y = l 0 ,z ) 
correlating factors: 

a 0 ~A U 0 (oj) 


kmM —y 7 

Fourier transformed object: 

F (k D - k m so) -A -J^O(K x , 7 -k, IQ 

(sign and factor k 2 defined at [3], eq. (12)) 


We then arrive at equation (46) in [3]. 

U£\x,y = l Q ,z ) = -tplloH JJ dK ^ Kz e^ l °d(K x , 7 


fc, lF z )e i[A ' x3:+A ' z2 ] ([3], equation (46)) 
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B. Changes to this Document 

[v2] March 2016 

- errata: 

- glossary entry <^(r): removed amplitude a as it is part of the complex phase 

- replaced i with j in text after equations 6.2, 6.13, and 6.23 

- updated references: 

- implementation of the backpropagation algorithm (updated [21]) 

- clarified missing-apple-core problem (added [23]) 

- removed journal abbreviations 

- replaced ’object function’ by ’scattering potential’ for more consistency with literature on 
scattering theory 

- deepened discussion on validity of Rytov approximation 

- changed equations for backpropagation algorithms 6.13 and 6.23 to include background cor¬ 
rection 

[v3] October 2016 

- errata: 

- glossary entry <p(r): added logarithm of amplitude a to complex phase 

- added missing negation in the discussion of the Rytov approximation 

- corrected equation 4.1: The inhomogeneous Helmhotz equation contains the scattering 
potential and not the delta distribution. 

- replaced equation 4.2: In 2D, the Green’s function is not oc exp(r)/r; replaced with the 
definition of the Green’s function. 
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